Inferência Bayesiana

Listas

Author

Micael Egídio Papa da Silva / 211029236

Published

24-06-2024

Lista 2

Questão 1

O seu professor chega na sala de aula e mostra uma moeda. Você suspeita que a moeda possa ser falsa e ter duas caras. Considere a priori probabilidades iguais para os eventos da moeda ser falsa ou ser honesta (i.e. uma moeda bem equilibrada).

  1. Calcule a sua probabilidade de obter cara num lançamento dessa moeda.

Para calcular a probabilidade de obter “cara” em um lançamento dessa moeda, considerando que a moeda pode ser honesta ou falsa. Suponha que:

  • \(H\) seja o evento “a moeda é honesta”.
  • \(F\) seja o evento “a moeda é falsa”.
  • \(C_{i}\) seja o evento “obter cara no i-ésimo lançamento”.

Assumindo pelo enunciado:

  • A priori, as probabilidades de a moeda ser honesta ou falsa são iguais: $P(H) = P(F) = 0.5 $.
  • Se a moeda é honesta, a probabilidade de obter cara é \(P(C_{1}|H) = 0.5\).
  • Se a moeda é falsa, a probabilidade de obter cara é \(P(C_{1}|F) = 1\) (pois tem duas caras).

A probabilidade de obter cara será dada por:

\(P(C_{1}) = P(C_{1}|H)P(H) + P(C_{1}|F)P(F)\)

Substituindo os valores fornecidos:

\(P(C_{1}) = (0.5 \times 0.5) + (1 \times 0.5)\)

\(P(C_{1}) = 0.25 + 0.5\)

\(P(C_{1}) = 0.75\)

Portanto, a probabilidade de obter cara em um lançamento dessa moeda é \(0.75\) ou 75%.

  1. Se o professor lançar a moeda e o resultado for cara, qual é agora a probabilidade dela ser falsa?

Queremos encontrar \(P(F|C_{1})\), a probabilidade de a moeda ser falsa dado que o resultado foi cara. Segue via Teorema de Bayes:

\(P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})}\)

Logo:

\(P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})} = \frac{1 \times 0.5}{0.75} = \frac{0.5}{0.75} = \frac{2}{3}\)

Portanto, a probabilidade de a moeda ser falsa, dado que o resultado do lançamento foi cara, é \(\frac{2}{3}\) ou aproximadamente 66.67%.

  1. Se o professor lançar a moeda n vezes e obter n caras, qual é a probabilidade dela ser falsa? Estude o comportamento desta probabilidade para n grande.

Assumimos:

  • \(C^n\) o evento “obter cara em todos os \(n\) lançamentos”.

Queremos calcular \(P(F|C^n)\), a probabilidade de a moeda ser falsa dado que todos os \(n\) lançamentos resultaram em cara.

Segue via Teorema de Bayes:

\(P(F|C^n) = \frac{P(C^n|F)P(F)}{P(C^n)}\)

Onde:

  • \(P(F) = 0.5\) (a priori, a moeda é falsa).
  • \(P(H) = 0.5\) (a priori, a moeda é honesta).
  • \(P(C^n|F) = 1\) (se a moeda é falsa, a probabilidade de obter cara em todos os lançamentos é 1).
  • \(P(C^n|H) = 0.5^n\) (se a moeda é honesta, a probabilidade de obter cara em todos os \(n\) lançamentos é \(0.5^n\)).

Para encontrar \(P(C^n)\), recorremos a probabilidade total:

\(P(C^n) = P(C^n|H)P(H) + P(C^n|F)P(F)\)

Substituindo os valores fornecidos:

\(P(C^n) = (0.5^n \times 0.5) + (1 \times 0.5)\)

\(P(C^n) = 0.5^{n+1} + 0.5\)

Agora, substituímos na fórmula de Bayes:

\(P(F|C^n) = \frac{1 \times 0.5}{0.5^{n+1} + 0.5}\)

$P(F|C^n) = # Inferência Bayesiana {.unnumbered}

Observa-se que valores grandes de \(n\), \(0.5^n\) se aproxima de 0 muito rapidamente, porque \(0.5^n\) decresce exponencialmente.

Portanto, para \(n\) grande, \(P(F|C^n)\) se aproxima de:

\(P(F|C^n) \approx \frac{1}{0 + 1} = 1\)

Isso significa que, conforme \(n\) aumenta, a probabilidade de a moeda ser falsa se aproxima de 1 (ou 100%).

Portanto, se o professor lançar a moeda \(n\) vezes e obter \(n\) caras, a probabilidade de a moeda ser falsa tende a 1 conforme \(n\) se torna grande.

  1. Se o professor lançar a moeda uma vez e obter cara, qual é a probabilidade do pr´oximo lançamento ser cara?

Após um lançamento resultando em “cara”, a probabilidade de a moeda ser falsa, \(P(F|C)\), foi calculada como \(\frac{2}{3}\). Consequentemente, a probabilidade de a moeda ser honesta, \(P(H|C)\), será:

\(P(H|C) = 1 - P(F|C) = 1 - \frac{2}{3} = \frac{1}{3}\)

Para a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento foi “cara” basta obter a soma das probabilidades condicionais de obter “cara” no próximo lançamento, dadas as duas possibilidades sobre a moeda (falsa ou honesta):

\(P(C_2|C_1) = P(C_2|F)P(F|C_1) + P(C_2|H)P(H|C_1)\)

Onde:

  • \(P(C_2|C_1)\) é a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento foi “cara”.
  • \(P(C_2|F) = 1\) (se a moeda é falsa, a probabilidade de obter “cara” no próximo lançamento é 1).
  • \(P(C_2|H) = 0.5\) (se a moeda é honesta, a probabilidade de obter “cara” no próximo lançamento é 0.5).
  • \(P(F|C_1) = \frac{2}{3}\) (a probabilidade de a moeda ser falsa após o primeiro “cara”).
  • \(P(H|C_1) = \frac{1}{3}\) (a probabilidade de a moeda ser honesta após o primeiro “cara”).

Substituindo os valores, temos:

\(P(C_2|C_1) = (1 \times \frac{2}{3}) + (0.5 \times \frac{1}{3})\)

\(P(C_2|C_1) = \frac{2}{3} + \frac{1}{6}\)

\(P(C_2|C_1) = \frac{4}{6} + \frac{1}{6}\)

\(P(C_2|C_1) = \frac{5}{6}\)

Portanto, a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento resultou em “cara”, é \(\frac{5}{6}\) ou aproximadamente 83.33%.

  1. Explique porque é falso neste contexto a afirmação ”os dois lançamentos da moeda são independentes”, e explique qual seria a afirmação correta.

No contexto descrito, a afirmação “os dois lançamentos da moeda são independentes” é falsa porque o resultado de cada lançamento influencia nossa crença sobre a natureza da moeda (se ela é honesta ou falsa). Essa influência altera a probabilidade dos resultados subsequentes.

Vale ressaltar que dois eventos \(A\) e \(B\) são independentes se e somente se:

\(P(A \cap B) = P(A) \cdot P(B)\)

Ou, de forma equivalente:

\(P(B|A) = P(B)\)

Se os lançamentos fossem verdadeiramente independentes, o resultado de um lançamento não alteraria a probabilidade do resultado do próximo lançamento. No entanto, neste contexto, o resultado do primeiro lançamento afeta nossa crença sobre a natureza da moeda, que por sua vez afeta a probabilidade do segundo lançamento. Ou seja , os lançamentos são dependentes porque \(P(C_2 | C_1) \neq P(C_2)\). Aqui, \(C_1\) e \(C_2\) representam o primeiro e o segundo lançamento resultando em cara, respectivamente.

A afirmação correta seria: “Os dois lançamentos da moeda são condicionalmente dependentes”.

Em resumo, o resultado de um lançamento influencia a probabilidade do próximo lançamento devido à incerteza sobre a natureza da moeda. Isso demonstra que os lançamentos não são independentes; são, na verdade, condicionalmente dependentes.

Questão 2

Seja \(y_{1} , y_{2} ,\dots, y_{n}\) uma amostra da distribuição de Bernoulli com probabilidade de sucesso θ e considere uma distribuição a priori uniforme para θ.

  1. Ache a distribuição a posteriori de θ e a sua média e variância.

Para uma amostra de Bernoulli, a verossimilhança é dada por:

\(P(y_1, y_2, \ldots, y_n | \theta) = \theta^{\sum_{i=1}^n y_i} (1 - \theta)^{n - \sum_{i=1}^n y_i}\)

Seja \(k = \sum_{i=1}^n y_i\), que é o número de sucessos na amostra. Assim, a função de verossimilhança pode ser escrita como:

\(P(y_1, y_2, \ldots, y_n | \theta) = \theta^k (1 - \theta)^{n - k}\)

A distribuição a priori uniforme para \(\theta\) no intervalo \([0, 1]\) é dada por:

\(P(\theta) = 1 \quad \text{para} \quad \theta \in [0, 1]\)

Distribuição a posteriori:

\(P(\theta | y_1, y_2, \ldots, y_n) \propto P(y_1, y_2, \ldots, y_n | \theta) \cdot P(\theta)\)

Substituindo as funções:

\(P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k} \cdot 1\)

\(P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k}\)

Essa é a forma da distribuição Beta. Portanto, a distribuição a posteriori de \(\theta\) é uma distribuição Beta com parâmetros \(\alpha = k + 1\) e \(\beta = n - k + 1\):

\(\theta | y_1, y_2, \ldots, y_n \sim \text{Beta}(k + 1, n - k + 1)\)

Para uma distribuição \(\text{Beta}(\alpha, \beta)\), a média e a variância são dadas por:

\(\mathbb{E}[\theta] = \frac{\alpha}{\alpha + \beta}\)

\(\text{Var}(\theta) = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}\)

Substituindo os parâmetros \(\alpha = k + 1\) e \(\beta = n - k + 1\):

Média:

\(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}\)

Variância:

\(\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{((k + 1) + (n - k + 1))^2 ((k + 1) + (n - k + 1) + 1)}\)

\(\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{(n + 2)^2 (n + 3)}\)

  1. Mostre que é possı́vel expressar a esperança a posteriori de θ da forma \((1−w)\times E(θ)+w\times \hat{θ̂}\), onde E(θ) e θ̂ são respectivamente a esperança a priori e a estimativa máximo verossı́mil de θ, e interprete este resultado.

A distribuição a posteriori de \(\theta\) foi encontrada como \(\text{Beta}(k + 1, n - k + 1)\), tendo \(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}\).

A priori, temos uma distribuição uniforme para \(\theta\) no intervalo \([0, 1]\), cuja esperança é:

\(\mathbb{E}[\theta] = \frac{1}{2}\)

Sabendo que a estimativa de máxima verossimilhança (MLE) de \(\theta\) para uma amostra Bernoulli é a proporção de sucessos na amostra:

\(\hat{\theta} = \frac{k}{n}\)

Basta substituir os valores encontrados:

\(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}\)

\(\mathbb{E}[\theta] = \frac{1}{2}\)

\(\hat{\theta} = \frac{k}{n}\)

Vamos expressar \(\frac{k + 1}{n + 2}\) na forma desejada:

\(\frac{k + 1}{n + 2} = \left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}\)

Para verificar, expandimos a expressão:

\(\left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}\)

\(= \frac{1}{2} - \frac{1}{n + 2} + \frac{2k}{n(n + 2)}\)

\(\frac{k + 1}{n + 2} = \frac{n}{2(n + 2)} + \frac{2k}{n(n + 2)}\)

Efetuando mais uma manipulação algébrica:

\(\frac{k + 1}{n + 2} = \frac{1}{2} \left( \frac{n}{n + 2} \right) + \left( \frac{2}{n + 2} \right) \times \frac{k}{n}\)

Logo, \(w = \frac{2}{n + 2}\) e \((1 - w) = \frac{n}{n + 2}\). Isso significa que a esperança a posteriori de \(\theta\) é uma média ponderada entre a esperança a priori e a estimativa de máxima verossimilhança, com os pesos dependendo do tamanho da amostra \(n\).

Ou seja:

  • Quando a amostra é pequena (\(n\) pequeno), a esperança a priori tem um peso maior na esperança a posteriori.

  • Quando a amostra é grande (\(n\) grande), a estimativa de máxima verossimilhança tem um peso maior, e a esperança a posteriori se aproxima da MLE.

Isso reflete a intuição de que com mais dados, nossa crença sobre \(\theta\) é mais fortemente influenciada pela evidência observada do que pelas crenças a priori.

  1. Se \(y_{n+1}\) é uma observação futura deste processo de Bernoulli, ache a distribuição preditiva \(p(y_{n+1} |y_{1} , y_{2} ,\dots, y_{n})\).

Distribuição preditiva de \(y_{n+1}\):

Para obter a distribuição preditiva de \(y_{n+1}\), precisamos calcular a probabilidade de \(y_{n+1}\) condicionado aos dados observados \(y_1, y_2, \dots, y_n\). A distribuição de \(y_{n+1}\) dado \(\theta\) é Bernoulli com parâmetro \(\theta\):

\(p(y_{n+1} | \theta) = \theta^{y_{n+1}} (1 - \theta)^{1 - y_{n+1}}\)

A distribuição preditiva é obtida integrando \(\theta\) fora desta expressão, ponderando pela distribuição a posteriori de \(\theta\), encontrada previamente:

\(p(y_{n+1} | y_1, y_2, \dots, y_n) = \int_0^1 p(y_{n+1} | \theta) p(\theta | y_1, y_2, \dots, y_n) \, d\theta\)

Substituindo \(p(y_{n+1} | \theta)\) e \(p(\theta | y_1, y_2, \dots, y_n)\):

\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)

\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)

Para \(y_{n+1} = 1\):

\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)

\(= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^{k + 1 - 1} (1 - \theta)^{n - k} \, d\theta\)

\(= \frac{B(k + 2, n - k + 1)}{B(k + 1, n - k + 1)}\)

\(= \frac{\frac{\Gamma(k + 2) \Gamma(n - k + 1)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}\)

\(= \frac{(k + 1) \Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 3)}{\Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 2)}\)

\(= \frac{(k + 1)}{n + 2}\)

Para \(y_{n+1} = 0\):

\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)

\(= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^k (1 - \theta)^{n - k + 1 - 1} \, d\theta\)

\(= \frac{B(k + 1, n - k + 2)}{B(k + 1, n - k + 1)}\)

\(= \frac{\frac{\Gamma(k + 1) \Gamma(n - k + 2)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}\)

\(= \frac{\Gamma(n - k + 2) / \Gamma(n + 3)}{\Gamma(n - k + 1) / \Gamma(n + 2)}\)

\(= \frac{(n - k + 1)}{n + 2}\)

Assim, a distribuição preditiva \(p(y_{n+1} | y_1, y_2, \dots, y_n)\) é:

\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \frac{k + 1}{n + 2}\)

\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \frac{n - k + 1}{n + 2}\)

A distribuição preditiva de uma observação futura \(y_{n+1}\) é a média ponderada das proporções de sucessos e fracassos observadas, ajustadas pelo tamanho da amostra, refletindo o princípio de suavização de Laplace, onde mesmo com uma grande quantidade de dados, há um pequeno ajuste baseado na crença a priori.

Questão 4

No exercício 2, calcule:

  1. a estimativa bayesiana para Perda Quadrática

Para calcular a estimativa Bayesiana sob a perda quadrática, precisamos encontrar o estimador que minimiza o valor esperado da perda quadrática, geralmente expressa como \(L(\theta, d) = (\theta - d)^2\), onde \(\theta\) é o parâmetro verdadeiro e \(d\) é a decisão ou estimativa tomada. Logo:

\(\mathbb{E}[(\theta - d)^2 | y_1, y_2, \dots, y_n]\)

onde \(\theta\) segue a distribuição a posteriori dado os dados observados \(y_1, y_2, \dots, y_n\).

Dado que \(\theta\) tem distribuição a posteriori \(\text{Beta}(k + 1, n - k + 1)\), a estimativa Bayesiana que minimiza a perda quadrática é a média da distribuição a posteriori, afinal, a média de uma distribuição é o estimador de mínimo erro quadrático médio (MSE).

A média da distribuição \(\text{Beta}(\alpha, \beta)\) é \(\frac{\alpha}{\alpha + \beta}\). Para a nossa distribuição a posteriori \(\text{Beta}(k + 1, n - k + 1)\):

\(\mathbb{E}[\theta | y_1, y_2, \dots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}\)

Portanto, a estimativa Bayesiana de \(\theta\) para a perda quadrática, dada a amostra de dados \([y_1, y_2, \dots, y_n]\), é:

\(d^* = \frac{k + 1}{n + 2}\)

A estimativa Bayesiana \(d^*\) corresponde à média a posteriori do parâmetro \(\theta\), levando em conta tanto a informação a priori (que \(\theta\) é uniformemente distribuído entre 0 e 1) quanto a informação proveniente da amostra (\(k\) sucessos em \(n\) tentativas). O ajuste de \(+1\) no numerador e \(+2\) no denominador é uma forma de regularização que evita extremos baseados em amostras pequenas, fornecendo uma suavização estatística que protege contra conclusões extremas com dados limitados, refletindo um compromisso entre o conhecimento a priori sobre a distribuição de \(\theta\) e a evidência obtida da amostra, alinhando-se com os princípios fundamentais da inferência Bayesiana.

  1. o limite da estimativa bayesiana para Perda Zero-Um quando \(\epsilon→ 0\). No caso especial que \(n = 12\), \(s = \sum ^{12}_{i=1} y_{i} = 9\)

A função de perda zero-um é dada por \(L(\theta, d) = 1\) se \(\theta \neq d\) e \(L(\theta, d) = 0\) se \(\theta = d\). Essencialmente, essa perda penaliza decisões incorretas sem considerar quão “erradas” elas estão, ao contrário da perda quadrática que penaliza de acordo com o quadrado da diferença entre a verdadeira \(\theta\) e a estimativa \(d\).

Cabe ressaltar que na estimativa para perda zero-um, procuramos minimizar o erro de classificação. Em geral, o estimador que minimiza a esperança da perda zero-um sob uma distribuição de probabilidades é o valor de \(\theta\) que maximiza a distribuição a posteriori, ou seja, a moda da distribuição a posteriori.

No nosso caso especial, a distribuição a posteriori de \(\theta\) após observar os dados é a \(\text{Beta}(s + 1, n - s + 1)\), onde \(s = \sum_{i=1}^n y_i\). Para \(n = 12\) e \(s = 9\), temos:

\(\theta | y_1, y_2, \dots, y_{12} \sim \text{Beta}(9 + 1, 12 - 9 + 1) = \text{Beta}(10, 4)\)

A moda de uma distribuição \(\text{Beta}(\alpha, \beta)\), para \(\alpha > 1\) e \(\beta > 1\), é dada por:

\(\text{Mode}[\theta] = \frac{\alpha - 1}{\alpha + \beta - 2}\)

Para a nossa distribuição \(\text{Beta}(10, 4)\):

\(\text{Mode}[\theta] = \frac{10 - 1}{10 + 4 - 2} = \frac{9}{12} = 0.75\)

O limite de \(\epsilon \rightarrow 0\) na pergunta sugere considerar um caso extremo da função de perda ou dos parâmetros envolvidos. Neste contexto, normalmente analisamos o comportamento da estimativa à medida que o tamanho da amostra cresce ou como a estimativa se comporta sob condições extremas da amostra ou dos parâmetros, todavia, como já calculamos a moda, o que seria o resultado para a perda zero-um, não há dependência direta de \(\epsilon\) no cálculo, a menos que \(\epsilon\) esteja relacionado a algum tipo de regularização ou parâmetro de suavização não especificado explicitamente na questão.

  1. a estimativa bayesiana sob Perda Absoluta

Para encontrar a estimativa Bayesiana sob a perda absoluta, precisamos minimizar o valor esperado da perda absoluta. A função de perda absoluta é dada por \(L(\theta, d) = |\theta - d|\), onde \(\theta\) é o parâmetro verdadeiro e \(d\) é a decisão ou estimativa tomada.

O estimador que minimiza a perda absoluta é a mediana da distribuição a posteriori, afinal, a mediana minimiza a soma das distâncias absolutas em relação a ela mesma.

A distribuição a posteriori de \(\theta\) após observar os dados $ y_1, y_2, , y_n $ é \(\text{Beta}(k + 1, n - k + 1)\), onde \(k = \sum_{i=1}^n y_i\). No caso especial onde $ n = 12 $ e $ k = s = 9 $, temos:

\(\theta | y_1, y_2, \ldots, y_{12} \sim \text{Beta}(10, 4)\)

Para uma distribuição \(\text{Beta}(\alpha, \beta)\), a mediana não tem uma fórmula fechada simples, mas pode ser aproximada ou calculada numericamente.

Uma aproximação comum para a mediana de uma distribuição \(\text{Beta}(\alpha, \beta)\) é:

\(\text{Mediana} \approx \frac{\alpha - \frac{1}{3}}{\alpha + \beta - \frac{2}{3}}\)

Substituindo \(\alpha = 10\) e \(\beta = 4\):

\(\text{Mediana} \approx \frac{10 - \frac{1}{3}}{10 + 4 - \frac{2}{3}} = \frac{9.6667}{13.3333} \approx 0.725\)

Portanto, a estimativa Bayesiana de \(\theta\) sob a perda absoluta, dada a amostra específica, é aproximadamente 0.725.

Ao optar por uma aproximação via python, obtemos:

Show/hide code
# Parâmetros da distribuição Beta
a, b = 10, 4

median = round(beta.median(a, b), 3)



print(f"A mediana da distribuição Beta({a}, {b}): {median}")
A mediana da distribuição Beta(10, 4): 0.725

De modo a encontrar o mesmo resultado para a estimativa Bayesiana de \(\theta\) sob a perda absoluta.

  1. um intervalo HPD com nível 99%.

O intervalo de alta densidade posterior (HPD - Highest Posterior Density) é um intervalo que contém uma determinada percentagem da densidade de probabilidade e onde todas as probabilidades dentro do intervalo são maiores ou iguais às probabilidades fora do intervalo. Para calcular um intervalo HPD de 99% para a distribuição a posteriori de \(\theta\), onde \(\theta\) segue uma distribuição \(\text{Beta}(10, 4)\), utilizaremos o software python, pois a distribuição Beta não permite uma solução analítica simples para o intervalo HPD.

Para \(n = 12\) e \(s = 9\), a distribuição a posteriori de \(\theta\) é \(\text{Beta}(10, 4)\). A fórmula da função de densidade para a distribuição Beta é:

\(f(\theta; \alpha, \beta) = \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}\)

onde \(B(\alpha, \beta)\) é a função Beta, e para nosso caso específico, \(\alpha = 10\) e \(\beta = 4\).

Logo:

Show/hide code
# Parâmetros da Beta
a, b = 10, 4


def hpd_interval(a, b, level=0.99):
    
    lower_bound = (1 - level) / 2
    upper_bound = 1 - lower_bound
    
    lower_quantile = beta.ppf(lower_bound, a, b)
    upper_quantile = beta.ppf(upper_bound, a, b)
    
    return lower_quantile, upper_quantile


hpd_99 = hpd_interval(a, b, 0.99)

hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))

print(f"Intervalo HPD de 99% para Beta({a}, {b}): {hpd_99_formatted}")
Intervalo HPD de 99% para Beta(10, 4): (0.379, 0.943)

O intervalo resultante explicita a região onde a densidade da distribuição Beta é mais alta e contém 99% da probabilidade total.

Questão 8

Suponha que \((x_{1} , x_{2} , x_{3})\) dado \(p_{1} , p_{2} , p_{3}\) segue uma distribuição Multinomial com parâmetros \(n\) e \((p_{1} , p_{2} , p_{3})\) , onde \(p_{i} ≥ 0 \space \text{e} \space p_{1} + p_{2} + p_{3} = 1\), e que, a priori, \((p_{1} , p_{2} , p_{3})\) segue uma distribuição de Dirichlet com parâmetros \((α_{1} , α_{2} , α_{3} )\):

  1. Ache a distribuição a posteriori de \(p_{1} , p_{2} , p_{3}\) e as distribuições a posteriori marginais de \(p_{i}\) (i = 1, 2, 3).

A distribuição a posteriori de \((p_1, p_2, p_3)\) é uma distribuição de Dirichlet com parâmetros atualizados.

A distribuição Dirichlet a priori é:

\((p_1, p_2, p_3) \sim \text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)\)

A função de verossimilhança para a distribuição multinomial é:

\(P(x_1, x_2, x_3 | p_1, p_2, p_3) = \frac{n!}{x_1! x_2! x_3!} p_1^{x_1} p_2^{x_2} p_3^{x_3}\)

Afinal, a probabilidade de obter exatamente \(x_1\) ocorrências do resultado 1, \(x_2\) ocorrências do resultado 2, e \(x_3\) ocorrências do resultado 3 é dada pelo produto das probabilidades individuais elevadas ao número de ocorrências, multiplicado pelo número de permutações possíveis.

A distribuição a posteriori será:

\(P(p_1, p_2, p_3 | x_1, x_2, x_3) \propto p_1^{\alpha_1 - 1} p_2^{\alpha_2 - 1} p_3^{\alpha_3 - 1} \cdot p_1^{x_1} p_2^{x_2} p_3^{x_3}\)

\(\propto p_1^{\alpha_1 + x_1 - 1} p_2^{\alpha_2 + x_2 - 1} p_3^{\alpha_3 + x_3 - 1}\)

Portanto, a distribuição a posteriori é uma distribuição de Dirichlet com parâmetros \((\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\):

\((p_1, p_2, p_3) | (x_1, x_2, x_3) \sim \text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\)

As distribuições marginais de \(p_i\) para uma distribuição de Dirichlet são Beta. Para \(p_1\):

\(p_1 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_1 + x_1, \alpha_2 + x_2 + \alpha_3 + x_3)\)

Similarmente, para \(p_2\):

\(p_2 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_2 + x_2, \alpha_1 + x_1 + \alpha_3 + x_3)\)

E para \(p_3\):

\(p_3 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_3 + x_3, \alpha_1 + x_1 + \alpha_2 + x_2)\)

  1. Calcule as estimativas bayesianas de \(p_{i}\) e de \(p_{j}\)\(p_{i}\) sob Perda Quadrática (i, j = 1, 2, 3, i 6= j).

Estimativas Bayesianas de \(p_i\):

Para a distribuição Dirichlet \(\text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)\), a esperança de \(p_i\) é:

\(\mathbb{E}[p_i | x_1, x_2, x_3] = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}\)

Portanto, as estimativas Bayesianas de \(p_1\), \(p_2\), e \(p_3\) são:

\(\hat{p_1} = \frac{\alpha_1 + x_1}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\) \(\hat{p_2} = \frac{\alpha_2 + x_2}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\) \(\hat{p_3} = \frac{\alpha_3 + x_3}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\)

Estimativa de \(p_j - p_i\):

Para calcular a estimativa Bayesiana de \(p_j - p_i\) sob a perda quadrática, usaremos a linearidade da esperança. A estimativa é simplesmente a diferença das esperanças de \(p_j\) e \(p_i\):

\(\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \mathbb{E}[p_j | x_1, x_2, x_3] - \mathbb{E}[p_i | x_1, x_2, x_3]\)

Portanto:

\(\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \frac{\alpha_j + x_j}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3} - \frac{\alpha_i + x_i}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\)

Podemos assim, concatenar os resultados obtidos :

  • A distribuição a posteriori de \((p_1, p_2, p_3)\) é \(\text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\).

  • As distribuições marginais de \(p_i\) são \(\text{Beta}(\alpha_i + x_i, \sum_{j \neq i} (\alpha_j + x_j))\).

  • As estimativas Bayesianas de \(p_i\) são \(\hat{p_i} = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}\).

  • As estimativas Bayesianas de \(p_j - p_i\) são \(\hat{p_j} - \hat{p_i} = \frac{\alpha_j + x_j}{\sum_{k=1}^3 (\alpha_k + x_k)} - \frac{\alpha_i + x_i}{\sum_{k=1}^3 (\alpha_k + x_k)}\).

Questão 11

É conhecido que 25% dos pacientes de um certo grupo que sofrem de enxaqueca melhoram após duas horas de serem tratados com um placebo. Para verificar se uma droga nova é melhor que o placebo, \(n = 20\) pacientes foram tratados com o placebo e verificou-se que após duas horas \(s = 8\) deles relataram ter melhorado. Seja \(\theta\) a probabilidade de um paciente tratado com a droga nova melhorar após duas horas.

  1. Especifique a hipótese nula \(H_0\) e a alternativa \(H_1\);
  • Hipótese nula \(H_0\): \(\theta = 0.25\) (a droga nova é tão eficaz quanto o placebo)
  • Hipótese alternativa \(H_1\): \(\theta > 0.25\) (a droga nova é mais eficaz que o placebo)
  1. Usando a distribuição a priori “não informativa” \(\theta \sim \text{Uniforme}(0, 1)\), calcule as chances relativas a priori e a posteriori de \(H_1\) e o correspondente Fator de Bayes.
  • \(P(\theta > 0.25)\)
  • \(P(\theta \leq 0.25)\)

A distribuição a priori uniforme \(\theta \sim \text{Uniforme}(0, 1)\):

\(P(\theta \leq 0.25) = 0.25\)

\(P(\theta > 0.25) = 1 - 0.25 = 0.75\)

As chances relativas a priori:

\(O(H_1) = \frac{P(\theta > 0.25)}{P(\theta \leq 0.25)} = \frac{0.75}{0.25} = 3\)

Sabemos que a distribuição a posteriori é \(\theta | s \sim \text{Beta}(9, 13)\).

Probabilidade a Posteriori de $ > 0.25 $:

\(P(\theta > 0.25 | s) = 1 - F_{\text{Beta}}(0.25 | 9, 13)\)

Consequentemente:

\(P(\theta \leq 0.25 | s) = F_{\text{Beta}}(0.25 | 9, 13)\)

Via python:

Show/hide code
# Parâmetros da distribuição Beta
a, b = 9, 13

posterior_prob_H1 = 1 - beta.cdf(0.25, a, b)
posterior_prob_H0 = beta.cdf(0.25, a, b)

print(f"P(θ > 0.25 | s) = {posterior_prob_H1:.3f}")
print(f"P(θ ≤ 0.25 | s) = {posterior_prob_H0:.3f}")
P(θ > 0.25 | s) = 0.944
P(θ ≤ 0.25 | s) = 0.056

Chances Relativas a Posteriori:

\(O(H_1 | s) = \frac{P(\theta > 0.25 | s)}{P(\theta \leq 0.25 | s)} = \frac{0.944}{0.056} = 16.810\)

Fator de Bayes:

\(B_{10}(s) = \frac{O(H_1 | s)}{O(H_1)} = \frac{16.810}{3} = 5.603\)

  1. Seja \(d = 1\) a decisão de rejeitar \(H_0\) e \(d = 0\) a de não rejeitar. Considere a função de perda de Neyman para a qual é 5 vezes mais custoso rejeitar \(H_0\) quando ela é verdadeira do que não rejeitar quando ela é falsa. Calcule a decisão ótima a posteriori;

Definindo a função de perda de Neyman:

  • \(L(d = 1, \theta \in H_0) = 5\)
  • \(L(d = 0, \theta \not\in H_0) = 5\)
  • \(L(d = 1, \theta \not\in H_0) = 0\)
  • \(L(d = 0, \theta \in H_0) = 0\)

Considerando a função de perda de Neyman:

  • \(a_0 = L(d = 1, \theta \leq 0.25) = 5L(d = 0, \theta \leq 0.25) = 5a_1\)

A decisão ótima é rejeitar \(H_0\) se:

\(P(\theta > 0.25 | s) > \frac{a_0}{a_0 + a_1} = \frac{5a_1}{5a_1 + a_1} = \frac{5}{6} = 0.833\)

Como \(P(\theta > 0.25 | s) = 0.944\) é maior do que \(0.833\), a decisão ótima é rejeitar \(H_0\).

  1. É razoável chamar essa distribuição a priori de “não informativa” nesse problema? Se a sua resposta for negativa, sugira uma outra distribuição a priori e refaça os cálculos anteriores.

Uma distribuição uniforme pode ser vista como “não informativa” em um certo sentido, mas pode não ser a mais adequada. Para hipóteses específicas como \(H_0\) e \(H_1\), uma distribuição Beta poderia ser mais informativa.

A priori uniforme favorece \(H_1\) no sentido que \(P(\theta > 0.25) = 3P(\theta \leq 0.25)\). Uma forma de corrigir isso é usar uma distribuição Beta com \(\beta = 3\alpha\), o que dá lugar a uma priori imprópria alternativa:

\(\theta \sim \text{Beta}(\alpha, 3\alpha)\)

Para obter uma variância maximizada, precisamos definir \(\alpha \to 0\). Então:

\(P(\theta) \propto \frac{1}{\theta(1 - \theta)}\)

Com a priori imprópria \(P(\theta) \propto \frac{1}{\theta(1 - \theta)}\), a distribuição a posteriori é proporcional a:

\(P(\theta | s) \propto P(s | \theta) P(\theta) \propto \theta^s (1 - \theta)^{n - s} \cdot \frac{1}{\theta(1 - \theta)} = \theta^{s-1} (1 - \theta)^{n-s-1}\)

Isto é uma distribuição \(\text{Beta}(s, n-s)\).

Para \(n = 20\) e \(s = 8\):

\(\theta | s \sim \text{Beta}(8, 12)\)

Show/hide code
alpha_prior, beta_prior = 1, 3
n = 20
s = 8

a_post = alpha_prior + s
b_post = beta_prior + n - s


posterior_prob_H1_beta = 1 - beta.cdf(0.25, a_post, b_post)
print(f"P(θ > 0.25 | s) = {posterior_prob_H1_beta:.3f}")
P(θ > 0.25 | s) = 0.904

A decisão de rejeitar \(H_0\) é mantida tanto com a priori uniforme quanto com a priori Beta(1, 3).

Lista 3.1

Questão 1

Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição de Poisson com média \(\theta\), e considere a priori que \(\theta\) tem uma distribuição Gama com parâmetros \(\alpha\) e \(\beta\) (isto é, com média \(\alpha\beta^{-1}\) e variância \(\alpha\beta^{-2}\)).

  1. Ache a distribuição a posteriori de \(\theta\) e a sua média e variância.

Dada uma amostra \(x_1, x_2, \ldots, x_n\) da distribuição de Poisson com média \(\theta\), a verossimilhança é:

\(P(x_1, x_2, \ldots, x_n | \theta) = \prod_{i=1}^n \frac{\theta^{x_i} e^{-\theta}}{x_i!}\)

\(P(x_1, x_2, \ldots, x_n | \theta) = \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \prod_{i=1}^n \frac{1}{x_i!}\)

Considerando que \(\theta\) tem uma distribuição Gama a priori com parâmetros \(\alpha\) e \(\beta\):

\(P(\theta) \propto \theta^{\alpha - 1} e^{-\beta \theta}\)

A distribuição a posteriori é:

\(P(\theta | x_1, x_2, \ldots, x_n) \propto P(x_1, x_2, \ldots, x_n | \theta) P(\theta)\)

Substituindo as funções de verossimilhança e a priori:

\(P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \cdot \theta^{\alpha - 1} e^{-\beta \theta}\)

\(P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{(\alpha + \sum_{i=1}^n x_i) - 1} e^{-(\beta + n) \theta}\)

Portanto, a distribuição a posteriori é uma distribuição Gama com parâmetros:

\(\alpha^* = \alpha + \sum_{i=1}^n x_i\)

\(\beta^* = \beta + n\)

A média e variância da distribuição a posteriori são:

\(\mathbb{E}[\theta | x] = \frac{\alpha^*}{\beta^*} = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}\)

\(\text{Var}[\theta | x] = \frac{\alpha^*}{(\beta^*)^2} = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}\)

  1. Mostre que é possível expressar a esperança a posteriori de \(\theta\) da forma \(w\bar{x} + (1-w)\alpha\beta^{-1}\), e interprete este resultado.

Sabemos que:

\(\mathbb{E}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}\)

Podemos escrever:

\(\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{\sum_{i=1}^n x_i}{\beta + n}\)

Definimos \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\), assim:

\(\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{n\bar{x}}{\beta + n}\)

Se definirmos \(w = \frac{n}{\beta + n}\), temos:

\(\mathbb{E}[\theta | x] = (1 - w) \frac{\alpha}{\beta} + w\bar{x}\)

Portanto, podemos escrever:

\(\mathbb{E}[\theta | x] = w\bar{x} + (1-w)\alpha\beta^{-1}\)

Isso mostra que a estimativa a posteriori é uma combinação ponderada entre a média a priori \(\alpha\beta^{-1}\) e a média amostral \(\bar{x}\).

  1. O que acontece na parte (ii) quando \(\beta\) é grande com \(\alpha\beta^{-1}\) fixo? Interprete!

Quando \(\beta\) é grande com \(\alpha\beta^{-1}\) fixo:

\(w = \frac{n}{\beta + n} \approx \frac{n}{\beta}\)

Se \(\beta \to \infty\), então \(w \to 0\), o que significa que a esperança a posteriori \(\mathbb{E}[\theta | x]\) será dominada pela média a priori \(\alpha\beta^{-1}\), resultando em:

\(\mathbb{E}[\theta | x] \approx (1 - 0) \alpha\beta^{-1} = \alpha\beta^{-1}\)

  1. Mostre que existe um número \(c\) tal que a variância a posteriori é maior do que a variância a priori sempre que \(\bar{x} > c\). Ache \(c\) e interprete este resultado.

A variância a posteriori é:

\(\text{Var}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}\)

A variância a priori é:

\(\text{Var}[\theta] = \frac{\alpha}{\beta^2}\)

Queremos encontrar \(c\) tal que \(\text{Var}[\theta | x] > \text{Var}[\theta]\) quando \(\bar{x} > c\).

\(\frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2} > \frac{\alpha}{\beta^2}\)

Multiplicando ambos os lados por \((\beta + n)^2 \beta^2\):

\((\alpha + \sum_{i=1}^n x_i) \beta^2 > \alpha (\beta + n)^2\)

\(\alpha \beta^2 + \sum_{i=1}^n x_i \beta^2 > \alpha (\beta^2 + 2\beta n + n^2)\)

\(\sum_{i=1}^n x_i \beta^2 > \alpha (2\beta n + n^2)\)

Dividindo ambos os lados por \(n \beta^2\):

\(\bar{x} > \frac{\alpha (2\beta n + n^2)}{n \beta^2}\)

\(\bar{x} > \frac{\alpha (2\beta + n)}{\beta^2}\)

Portanto, o valor crítico \(c\) é:

\(c = \frac{\alpha (2\beta + n)}{\beta^2} = \frac{\alpha}{\beta} \times (2 + \frac{n}{\beta}) = (2 + \frac{n}{\beta}) \times \mathbb{E}[\theta]\)

Questão 3

  1. Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição de Poisson com média \(\theta\), e considere a priori que \(\theta\) tem uma distribuição Gama com parâmetros \(\alpha = 1\) e \(\beta = 1\). Ache:
  1. A estimativa bayesiana de \(\theta\) no caso de perda quadrática.

A estimativa sob perda quadrática é imediata: \(\mathbb{E}[\theta | D] = \alpha^* / \beta^* = (s + 1) / (n + 1)\).

  1. O limite da estimativa bayesiana sob perda zero-um quando \(\epsilon \to 0\).

O limite da estimativa sob perda zero-um é a moda da distribuição a posteriori, que nesse caso é \((\alpha^* - 1) / \beta^* = s / (n + 1)\).

Agora, quando \(n = 10\) e \(\bar{x} = 1.55\), temos que \(\alpha^* = 16.50\) e \(\beta^* = 11\), de forma que a estimativa sob PA é a mediana \(\text{med}(\theta | D) \approx 1.470\).

Para o caso \(n = 10\) e \(\bar{x} = 1.55\), ache:

  1. A estimativa bayesiana sob perda absoluta.

Temos que \(\alpha^* = 16.50\) e \(\beta^* = 11\), de forma que a estimativa sob PA é a mediana \(\text{med}(\theta | D) \approx 1.470\).

  1. O intervalo HPD para \(\theta\) com nível 95%.

Vamos implementar os cálculos em Python para verificar os resultados obtidos anteriormente e, por fim, explicitar o intervalo HPD:

Show/hide code
n = 10
x_bar = 1.55
s = n * x_bar

alpha_post = 1 + s
beta_post = 1 + n

theta_quadratic = alpha_post / beta_post
print(f"Estimativa bayesiana sob perda quadrática: {theta_quadratic:.3f}")

theta_zero_one = (alpha_post - 1) / beta_post
print(f"Limite da estimativa sob perda zero-um: {theta_zero_one:.3f}")

theta_absolute = gamma.median(alpha_post, scale=1/beta_post)
print(f"Estimativa bayesiana sob perda absoluta: {theta_absolute:.3f}")

hpd_interval = gamma.interval(0.95, alpha_post, scale=1/beta_post)
print(f"Intervalo HPD de 95%: ({hpd_interval[0]:.3f}, {hpd_interval[1]:.3f})")
Estimativa bayesiana sob perda quadrática: 1.500
Limite da estimativa sob perda zero-um: 1.409
Estimativa bayesiana sob perda absoluta: 1.470
Intervalo HPD de 95%: (0.866, 2.306)

Questão 5

Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\) conhecida, e considere a distribuição a priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\).

  1. Ache a distribuição a posteriori de \(\mu\).

Dada a amostra \(x_1, x_2, \ldots, x_n\) da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\), a verossimilhança é:

\(f(x_i | \mu) = \sqrt{\frac{\phi}{2\pi}} \exp\left( -\frac{\phi}{2}(x_i - \mu)^2 \right)\)

A função de verossimilhança conjunta para a amostra é:

\(L(\mu | x_1, x_2, \ldots, x_n) = \left(\frac{\phi}{2\pi}\right)^{n/2} \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right)\)

A priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\) é dada por:

\(f(\mu) = \sqrt{\frac{\tau}{2\pi}} \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)\)

Onde teremos a distribuição a posteriori:

\(f(\mu | x_1, x_2, \ldots, x_n) \propto L(\mu | x_1, x_2, \ldots, x_n) f(\mu)\)

Substituindo as expressões:

\(f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right) \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)\)

Agrupando os termos quadráticos em \(\mu\):

\(\propto \exp\left( -\frac{1}{2} \left[ \phi \sum_{i=1}^n x_i^2 + \tau \mu_0^2 - 2\mu \left( \phi \sum_{i=1}^n x_i + \tau \mu_0 \right) + \mu^2 (\phi n + \tau) \right] \right)\)

Completando o quadrado:

\(f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{1}{2} \left[ (\phi n + \tau)\left(\mu - \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\right)^2 \right] \right)\)

Portanto, a distribuição a posteriori é uma Normal com os parâmetros atualizados:

\(\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)\)

  1. Mostre que é possível expressar a esperança a posteriori de \(\mu\) da forma \(w\bar{x} + (1 - w)\mu_0\), e interprete este resultado.

A esperança a posteriori é:

\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)

Podemos reescrever essa expressão como uma combinação ponderada entre a média amostral \(\bar{x}\) e a média a priori \(\mu_0\). Definimos \(w = \frac{\phi n}{\phi n + \tau}\):

\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = w \bar{x} + (1 - w) \mu_0\)

onde \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\).

Este resultado mostra que a estimativa a posteriori é uma média ponderada entre a média da amostra \(\bar{x}\) e a média a priori \(\mu_0\), com pesos \(w\) e \(1 - w\) respectivamente. O peso \(w\) depende do tamanho da amostra e dos parâmetros \(\phi\) e \(\tau\).

  1. Se \(\bar{x}_m\) é a média de \(m\) observações futuras \(x_{n+1}, \ldots, x_{n+m}\), condicionalmente independentes de \(x_1, \ldots, x_n\), ache a distribuição preditiva \(p(\bar{x}_m | x_1, \ldots, x_n)\).

Se \(\bar{x}_m\) é a média de \(m\) observações futuras \(x_{n+1}, \ldots, x_{n+m}\), condicionalmente independentes de \(x_1, \ldots, x_n\), a distribuição preditiva é:

\(p(\bar{x}_m | x_1, \ldots, x_n)\)

A média de \(m\) observações futuras é normalmente distribuída com a mesma média \(\mu\) e variância \(\frac{1}{\phi m}\). Portanto, a média preditiva é uma média ponderada da distribuição posterior de \(\mu\) e da distribuição das observações futuras.

A variância total é a soma das variâncias individuais:

\(\text{Var}(\bar{x}_m | x_1, \ldots, x_n) = \text{Var}(\mu | x_1, \ldots, x_n) + \text{Var}(x_{n+1}, \ldots, x_{n+m} | \mu)\)

\(= \frac{1}{\phi n + \tau} + \frac{1}{\phi m}\)

A distribuição preditiva de \(\bar{x}_m\) é:

\(\bar{x}_m | x_1, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi n \bar{x} + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} + \frac{1}{\phi m} \right)\)

  1. Discuta o que acontece com os resultados anteriores quando a distribuição a priori \(p(\mu) \propto 1\) (i.e., o caso limite quando \(\tau \to 0\)).

Quando \(\tau \to 0\), a priori se torna não informativa, ou seja, \(p(\mu) \propto 1\). Nesse caso, o peso \(w\) na expressão da esperança a posteriori se torna 1, pois:

\(w = \frac{\phi n}{\phi n + \tau} \approx 1\)

Portanto, a esperança a posteriori se torna a média amostral:

\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] \approx \bar{x}\)

E a variância a posteriori se torna:

\(\text{Var}(\mu | x_1, x_2, \ldots, x_n) \approx \frac{1}{\phi n}\)

Isso mostra que, com uma priori não informativa, a distribuição a posteriori é dominada pelos dados da amostra, e a média a posteriori é simplesmente a média amostral.

Questão 6

Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\) conhecida, e considere a distribuição a priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\).

(a) Ache o estimador bayesiano de \(\mu\) no caso de (i) perda quadrática \((L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2)\), (ii) perda absoluta \((L(\hat{\mu}, \mu) = |\hat{\mu} - \mu|)\) e (iii) perda zero-um \((L(\hat{\mu}, \mu) = I(|\hat{\mu} - \mu| \geq \epsilon))\).

Para uma função de perda quadrática \(L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2\), o estimador bayesiano é a média da distribuição a posteriori.

Como sabemos da questão anterior que a distribuição a posteriori é uma Normal com parâmetros atualizados:

\(\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)\)

A média dessa distribuição é:

\(\hat{\mu}_{\text{PQ}} = \mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)

Para a distribuição normal, a média, a mediana e a moda são todos iguais devido à sua simetria. Assim, os estimadores sob perda quadrática (i), perda absoluta (ii) e perda zero-um (iii) são todos iguais a:

\(\hat{\mu} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)

Essa equivalência é uma consequência direta da simetria da distribuição normal, onde a distribuição é perfeitamente simétrica em torno da média, que coincide com a mediana e a moda.

(b) Ache o intervalo HPD para \(\mu\) com nível \(100(1 - \alpha)\%\).

O intervalo HPD (Highest Posterior Density) para uma distribuição Normal \(\mathcal{N}(\mu_{\text{post}}, \sigma^2_{\text{post}})\) é encontrado nos quantis da distribuição que correspondem ao nível de confiança desejado.

Os limites do intervalo HPD de \(100(1 - \alpha)\%\) são:

\(\left[ \mu_{\text{post}} - z_{\alpha/2} \sigma_{\text{post}}, \mu_{\text{post}} + z_{\alpha/2} \sigma_{\text{post}} \right]\)

onde \(\mu_{\text{post}} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\) e \(\sigma^2_{\text{post}} = \frac{1}{\phi n + \tau}\).

Vamos calcular esses valores:

Show/hide code
n = 10
x_bar = 1.55
phi = 1
tau = 1
mu_0 = 0

sum_x = n * x_bar
mu_post = (phi * sum_x + tau * mu_0) / (phi * n + tau)
sigma_post = np.sqrt(1 / (phi * n + tau))

alpha = 0.05
z = norm.ppf(1 - alpha / 2)
hpd_interval = [mu_post - z * sigma_post, mu_post + z * sigma_post]

print(f"Intervalo HPD de 95%: [{hpd_interval[0]:.3f}, {hpd_interval[1]:.3f}]")
Intervalo HPD de 95%: [0.818, 2.000]

Lista 3.2

Questão 1

Considere o modelo \(X \sim \text{Uniforme} (\theta, \theta + 1)\) \((-\infty < \theta < \infty)\).

(a) Mostre que \(\theta\) é um parâmetro de locação.

Um parâmetro \(\theta\) é um parâmetro de locação se, ao adicionar uma constante \(c\) a \(\theta\), a distribuição resultante for uma translação da distribuição original.

Para o modelo \(X \sim \text{Uniforme} (\theta, \theta + 1)\):

\[ f_X(x; \theta) = \begin{cases} 1 & \text{se } \theta \leq x \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Se adicionarmos uma constante \(c\) a \(\theta\), a nova distribuição será \(X \sim \text{Uniforme} (\theta + c, \theta + c + 1)\):

\[ f_X(x; \theta + c) = \begin{cases} 1 & \text{se } \theta + c \leq x \leq \theta + c + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Essa é uma translação da distribuição original, mostrando que \(\theta\) é um parâmetro de locação.

(b) Dada a priori não informativa usual para o modelo de locação, \(p(\theta) \propto 1\), e uma amostra aleatória \(X_1, \ldots, X_n\) do modelo acima, discuta a propriedade da distribuição a posteriori.

A priori não informativa usual para um modelo de locação é \(p(\theta) \propto 1\).

Dada uma amostra aleatória \(X_1, \ldots, X_n\) do modelo \(\text{Uniforme} (\theta, \theta + 1)\):

A verossimilhança é:

\(L(\theta | X_1, \ldots, X_n) \propto \prod_{i=1}^n f_X(x_i; \theta)\)

\[ L(\theta | X_1, \ldots, X_n) = \begin{cases} 1 & \text{se } \theta \leq \min(X_i) \text{ e } \max(X_i) \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Dado que a priori é \(p(\theta) \propto 1\), a distribuição a posteriori é proporcional à verossimilhança:

\(p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n)\)

\[ p(\theta | X_1, \ldots, X_n) = \begin{cases} \frac{1}{1 - (\max(X_i) - \min(X_i))} & \text{se } \max(X_i) - 1 \leq \theta \leq \min(X_i) \\ 0 & \text{caso contrário} \end{cases} \]

Portanto, a distribuição a posteriori de \(\theta\) é uma distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\).

(c) Na situação da parte (b), ache:

  1. os estimadores bayesianos sob Perda Quadrática e sob Perda Absoluta;
  • Perda Quadrática: Para a função de perda quadrática \(L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2\), o estimador bayesiano é a média da distribuição a posteriori.

Dado que a distribuição a posteriori é uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\), a média é:

\(\hat{\theta}_{\text{PQ}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}\)

  • Perda Absoluta: Para a função de perda absoluta \(L(\hat{\theta}, \theta) = |\hat{\theta} - \theta|\), o estimador bayesiano é a mediana da distribuição a posteriori.

A mediana de uma distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\) é:

\(\hat{\theta}_{\text{PA}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}\)

  1. o intervalo HPD para \(\theta\) com credibilidade \(100(1 - \alpha)\%\);

Para a distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\), o intervalo HPD \(100(1 - \alpha)\%\) é simplesmente o intervalo central que cobre \(100(1 - \alpha)\%\) da distribuição:

\(\text{HPD} = \left[ X_n - 1 + \frac{\alpha}{2} (X_1 - (X_n - 1)), X_1 - \frac{\alpha}{2} (X_1 - (X_n - 1)) \right]\)

Onde \(\min(X_i) = X_1\) e \(\max(X_i) = X_n\)

  1. a probabilidade a posteriori do evento \(\theta \geq \theta_0\).

Para encontrar a probabilidade a posteriori do evento \(\theta \geq \theta_0\):

\(P(\theta \geq \theta_0 | X_1, \ldots, X_n)\)

Para a distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\):

Se \(\theta_0 < \max(X_i) - 1\), então $ P(_0 | X_1, , X_n) = 1$.

Se \(\theta_0 > \min(X_i)\), então $ P(_0 | X_1, , X_n) = 0$.

Se \(\max(X_i) - 1 \leq \theta_0 \leq \min(X_i)\):

\(P(\theta \geq \theta_0 | X_1, \ldots, X_n) = \frac{\min(X_i) - \theta_0}{\min(X_i) - (\max(X_i) - 1)}\)

Questão 2

Considere o modelo \(X \sim \text{Uniforme} (0, \theta)\) \((0 < \theta < \infty)\).

  1. Mostre que \(\sigma = \theta^{-1}\) é um parâmetro de escala.

Considere \(X \sim \text{Uniforme}(0, \theta)\). A função densidade de probabilidade (pdf) de \(X\) é:

\(f_X(x; \theta) = \theta^{-1} \mathbb{I}_{\{0 \leq x \leq \theta\}}\)

Vamos redefinir a distribuição em termos de \(\sigma\). Como \(\sigma = \theta^{-1}\), então \(\theta = \sigma^{-1}\). A pdf em termos de \(\sigma\) é:

\(f_X(x; \sigma) = \sigma \mathbb{I}_{\{0 \leq x \leq \sigma^{-1}\}}\)

Vamos multiplicar \(\sigma\) por uma constante \(c\). A nova densidade será:

\(f_Y(y; c\sigma) = c\sigma \mathbb{I}_{\{0 \leq y \leq (c\sigma)^{-1}\}}\)

Essa é uma escala da densidade original, mostrando que \(\sigma = \theta^{-1}\) é um parâmetro de escala.

  1. Dada a priori não informativa usual para o modelo de escala, \(p(\sigma) \propto \sigma^{-1}\), e uma amostra aleatória \(X_1, \ldots, X_n\) do modelo acima, discuta a propriedade da distribuição a posteriori.

Segue que:

\(L(\theta | X_1, \ldots, X_n) = \prod_{i=1}^n f_X(x_i; \theta)\)

Como \(X_i \sim \text{Uniforme}(0, \theta)\), a verossimilhança é:

\(L(\theta | X_1, \ldots, X_n) = \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}}\)

A priori não informativa é \(p(\sigma) \propto \sigma^{-1}\). Como \(\sigma = \theta^{-1}\), temos \(p(\theta) \propto \theta^{-1}\). Logo, a distribuição posteriori será:

\(p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n) p(\theta)\)

\(p(\theta | X_1, \ldots, X_n) \propto \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}} \theta^{-1}\)

\(p(\theta | X_1, \ldots, X_n) \propto \theta^{-(n+1)} \mathbb{I}_{\{\theta \geq \max(X_i)\}}\)

Por fim:

\(p(\theta | X_1, \ldots, X_n) = \frac{(n \max(X_i))^n}{\theta^{n+1}} \mathbb{I}_{\{\theta \geq \max(X_i)\}}\)

Esta é uma distribuição de Pareto com forma \(n\) e escala \(\max(X_i)\).

  1. Na situação da parte (b), ache:

  2. os estimadores bayesianos sob Perda Quadrática e sob Perda Absoluta.

  • Perda Quadrática: Para a função de perda quadrática \(L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2\), o estimador bayesiano é a média da distribuição a posteriori. A média de uma distribuição Pareto \(\text{Pareto}(\max(X_i), n)\) é dada por: \(\hat{\theta}_{\text{PQ}} = \frac{n \max(X_i)}{n - 1}\) , para \(n > 1\).

  • Perda Absoluta: Para a função de perda absoluta \(L(\hat{\theta}, \theta) = |\hat{\theta} - \theta|\), o estimador bayesiano é a mediana da distribuição a posteriori. A mediana de uma distribuição Pareto \(\text{Pareto}(\max(X_i), n)\) é dada por:

\(\hat{\theta}_{\text{PA}} = \max(X_i) \cdot 2^{1/(n)}\)

  1. o intervalo HPD para \(\theta\) com credibilidade \(100(1 - \alpha)\%\).

Para uma distribuição Pareto \(\text{Pareto}(\max(X_i), n+1)\), o intervalo HPD \(100(1 - \alpha)\%\) pode ser calculado encontrando os quantis apropriados. Seja \(F_{\text{Pareto}}^{-1}(p)\) a função quantil da distribuição Pareto:

\(\text{HPD} = [ \max(X_i), \max(X_i) (1 - \alpha)^{-\frac{1}{n+1}} ]\)

  1. a probabilidade a posteriori do evento \(\theta > \theta_0\).

A probabilidade a posteriori do evento \(\theta > \theta_0\) é dada pela função de sobrevivência da distribuição Pareto:

\(P(\theta > \theta_0 | X_1, \ldots, X_n) = \left( \frac{\max(X_i)}{\theta_0} \right)^{n} \mathbb{I}_{\{\theta_0 \geq \max(X_i)\}}\)

Questão 3

Considere o modelo \(x_1, \ldots, x_n | \theta \sim \text{Ber}(\theta)\).

  1. Calcule a priori de Jeffreys e mostre que ela é própria.

A priori de Jeffreys é uma escolha não informativa que é proporcional à raiz quadrada do determinante da matriz de informação de Fisher. Para o modelo Bernoulli \(\text{Ber}(\theta)\), vamos calcular a matriz de informação de Fisher e encontrar a priori de Jeffreys.

Para uma amostra \(x_1, x_2, \ldots, x_n\) de uma distribuição Bernoulli com parâmetro \(\theta\), a função de verossimilhança é:

\(L(\theta) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}\)

A log-verossimilhança é:

\(\log L(\theta) = \left(\sum_{i=1}^n x_i\right) \log \theta + \left(n - \sum_{i=1}^n x_i\right) \log (1 - \theta)\)

A informação de Fisher é a expectativa do quadrado da derivada da log-verossimilhança em relação a \(\theta\):

\(I(\theta) = -\mathbb{E}\left[\frac{\partial^2}{\partial \theta^2} \log L(\theta)\right]\)

Calculando a segunda derivada da log-verossimilhança:

\(\frac{\partial}{\partial \theta} \log L(\theta) = \frac{\sum_{i=1}^n x_i}{\theta} - \frac{n - \sum_{i=1}^n x_i}{1 - \theta}\)

\(\frac{\partial^2}{\partial \theta^2} \log L(\theta) = -\frac{\sum_{i=1}^n x_i}{\theta^2} - \frac{n - \sum_{i=1}^n x_i}{(1 - \theta)^2}\)

Tomando a expectativa, considerando que \(\mathbb{E}[x_i] = \theta\):

\(I(\theta) = n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)\)

A priori de Jeffreys é proporcional à raiz quadrada da informação de Fisher:

\(\pi_J(\theta) \propto \sqrt{I(\theta)}\)

\(\pi_J(\theta) \propto \sqrt{n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)}\)

\(\pi_J(\theta) \propto \sqrt{\frac{n}{\theta(1 - \theta)}}\)

\(\pi_J(\theta) \propto \frac{1}{\sqrt{\theta(1 - \theta)}}\)

Esta é a distribuição Beta \(\text{Beta}(1/2, 1/2)\), que é própria porque a integral da densidade de Jeffreys sobre o intervalo \([0, 1]\) é finita.

  1. Considere uma amostra aleatória \(X_1, \ldots, X_n\) do modelo acima e a priori de Jeffreys. Ache o estimador bayesiano de \(\theta\) sob Perda Quadrática e explique como achar um intervalo HPD para \(\theta\) com credibilidade \(100(1 - \alpha)\%\).

Temos que a verossimilhança é:

\(L(\theta | x_1, \ldots, x_n) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}\)

Com a priori de Jeffreys, a posteriori resulta em:

\(\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i} \cdot \frac{1}{\sqrt{\theta(1 - \theta)}}\)

\(\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i - 1/2} (1 - \theta)^{n - \sum_{i=1}^n x_i - 1/2}\)

Esta é uma distribuição Beta \(\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)\).

O estimador bayesiano sob perda quadrática é a média da distribuição a posteriori Beta:

\(\hat{\theta}_{\text{PQ}} = \mathbb{E}[\theta | x_1, \ldots, x_n] = \frac{\sum_{i=1}^n x_i + 1/2}{n + 1}\)

Para encontrar o intervalo HPD, precisamos calcular os quantis da distribuição Beta posteriori \(\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)\).

Como a questão não explicita o total de observações e a soma das observações, podemos visualizar o comportamento do HPD ao passo em que o tamanho amostral aumenta:

Show/hide code
def hpd_interval(a, b, level=0.99, sample_size=100000):
    sample = beta.rvs(a, b, size=sample_size)
    sample_sorted = np.sort(sample)
    interval_idx_inc = int(np.floor(level * sample_size))
    n_intervals = sample_size - interval_idx_inc
    interval_width = sample_sorted[interval_idx_inc:] - sample_sorted[:n_intervals]
    min_idx = np.argmin(interval_width)
    hpd_min = sample_sorted[min_idx]
    hpd_max = sample_sorted[min_idx + interval_idx_inc]
    return hpd_min, hpd_max

def update_plot(n):
    media_x = 0.7  # Supondo uma média dos xi
    sum_x_i = media_x * n
    a = sum_x_i + 0.5
    b = n - sum_x_i + 0.5
    
    x = np.linspace(0, 1, 1000)
    y = beta.pdf(x, a, b)
    
    hpd_99 = hpd_interval(a, b, level=0.99)
    hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))
    
    plt.figure(figsize=(6, 3))
    plt.plot(x, y, label=f'Beta({a:.2f}, {b:.2f})')
    plt.fill_between(x, y, where=(x >= hpd_99[0]) & (x <= hpd_99[1]), color='gray', alpha=0.5, label='HPD')
    plt.title(f'Distribuição Beta e Intervalo HPD de 99%\nIntervalo HPD: {hpd_99_formatted}')
    plt.xlabel('θ')
    plt.ylabel('Densidade de Probabilidade')
    plt.legend()
    plt.show()

n_slider = IntSlider(min=1, max=100, step=1, value=10, description='n')

interactive_plot = interactive(update_plot, n=n_slider)

display(interactive_plot)

O gráfico supracitado foi feito pressupondo que a média de \(x_i\) fosse 0.7 , pode-se perceber que a média da distribuição posteriori se torna mais estável e menos sensível a pequenas alterações nos dados conforme 𝑛 aumenta. Isso significa que a estimativa pontual de 𝜃 se torna mais confiável.

Os parâmetros \(a\) e \(b\) da distribuição Beta aumentam linearmente com \(n\), tornando a distribuição mais estreita. Para tamanhos amostrais pequenos, a distribuição Beta toma formas largas e até assimétricas, dependendo dos valores observados. Conforme \(n\) aumenta, a distribuição tende a se tornar mais simétrica e mais próxima de uma distribuição normal centrada na média posteriori e ,consequentemente, reduzindo o comprimento do intervalo HPD, convergindo para um valor fixo muito próximo da média verdadeira de \(\theta\).

Lista 4

Questão 1

Considere uma amostra \(y_1, \ldots, y_n\) da distribuição Normal com média \(\mu\) e variância \(\sigma^2 = 1/\tau\) desconhecidas, e suponha que a priori a distribuição de \((\mu, \tau)\) é especificada da seguinte forma: \(\mu | \tau \sim N(\mu_0, 1/\lambda_0 \tau)\) e \(\tau \sim \text{Gama}(\alpha_0, \beta_0)\), onde \(\lambda_0\), \(\alpha_0\) e \(\beta_0\) são positivas.

  1. Ache a distribuição a posteriori de \(p(\mu, \tau | D)\) e as distribuições marginais \(p(\mu | D)\) e \(p(\tau | D)\).

Ressaltando:

  • Priori: \(\mu | \tau \sim N(\mu_0, 1/\lambda_0 \tau)\) \(\tau \sim \text{Gama}(\alpha_0, \beta_0)\)

A função de verossimilhança para uma amostra \(y_1, \ldots, y_n\) é:

\(p(y | \mu, \tau) = \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right)\)

Temos a posteriori:

\(p(\mu, \tau | y) \propto p(y | \mu, \tau) p(\mu | \tau) p(\tau)\)

Substituindo as expressões:

\[ \begin{align*} p(\mu, \tau \mid y) &\propto \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right) \\ &\quad \times \left( \frac{\lambda_0 \tau}{2\pi} \right)^{1/2} \exp\left( -\frac{\lambda_0 \tau}{2} (\mu - \mu_0)^2 \right) \\ &\quad \times \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)} \tau^{\alpha_0 - 1} \exp(-\beta_0 \tau) \end{align*} \]

Agrupando os termos dependentes de \(\mu\) e \(\tau\):

\(p(\mu, \tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp\left( -\frac{\tau}{2} \left[ \sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2 \right] - \beta_0 \tau \right)\)

Para integrar \(\mu\), completamos o quadrado no expoente:

\(\sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2\)

Expandindo e agrupando os termos:

\(\sum_{i=1}^n y_i^2 - 2\mu \sum_{i=1}^n y_i + n\mu^2 + \lambda_0 \mu_0^2 - 2\lambda_0 \mu \mu_0 + \lambda_0 \mu^2\)

Agrupando os termos quadráticos, lineares e constantes em \(\mu\):

\(\left( n + \lambda_0 \right) \mu^2 - 2 \left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right) \mu + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2\)

Completando o quadrado:

\(\left( n + \lambda_0 \right) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \left( \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 \right) - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}\)

Assim, o expoente se torna:

\(-\frac{\tau}{2} \left[ (n + \lambda_0) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0} \right]\)

A integral em \(\mu\) da expressão acima é uma integral gaussiana, logo:

\(\int_{-\infty}^{\infty} \exp \left( -\frac{\tau (n + \lambda_0)}{2} \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 \right) d\mu = \sqrt{\frac{2\pi}{\tau (n + \lambda_0)}}\)

Portanto, ao desprezarmos essa constante, a distribuição marginal de \(\tau\) é:

\(p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}}{2} + \beta_0 \right) \right)\)

Simplificando a constante no expoente:

\(p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n (y_i - \bar{y})^2}{2} + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{2(\lambda_0 + n)} + \beta_0 \right) \right)\)

Reconhecendo a forma da distribuição Gama, temos:

\(\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)\)

Dada a \(\tau\), a distribuição condicional de \(\mu\) é:

\(\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)\)

\(\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)\)

Dada a \(\tau\), a distribuição condicional de \(\mu\) é:

\(\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)\)

  1. Discuta neste contexto o uso da distribuição a priori não informativa \(p(\mu, \tau) \propto 1/\tau\). A distribuição a posteriori é própria? Qual é a relação da distribuição \(p(\mu | D)\) e os resultados clássicos?

Essa priori é escolhida para ser não informativa, significando que temos pouca ou nenhuma informação a priori sobre \(\mu\) e \(\tau\).

  • Distribuição a posteriori: Com a priori não informativa, a posteriori ainda é própria, pois as integrais de normalização convergem.

  • Comparação com Resultados Clássicos:

    • A distribuição marginal de \(\mu\) se aproxima da distribuição normal clássica para a média de uma amostra, especialmente à medida que \(n\) aumenta.

    • A distribuição marginal de \(\tau\) se aproxima da distribuição qui-quadrado clássica (quando reescalada), que é usada para inferência sobre a variância na estatística frequentista.

Questão 2c

Seja \(n = (n_1, \ldots, n_k)\) um vetor aleatório com distribuição multinomial e densidade \(p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}\), onde \(\theta = (\theta_1, \ldots, \theta_k)\), \(\theta_i > 0\) e \(\sum_{i=1}^k \theta_i = 1\). Considere a priori para \(\theta\) uma distribuição de Dirichlet com parâmetro \(\alpha = (\alpha_1, \ldots, \alpha_k)\), isto é \(p(\theta) \propto \prod_{i=1}^k \theta_i^{\alpha_i - 1}\).

No caso particular \(\alpha = (0, 0, \ldots, 0)\) a distribuição a priori de \(\theta\) é imprópria. Mostre que a distribuição a posteriori é própria se e somente se \(n_i > 0\) para \(i = 1, 2, \ldots, k\).

Quando \(\alpha = (0, 0, \ldots, 0)\), a distribuição a priori de \(\theta\) é:

\(p(\theta) \propto \prod_{i=1}^k \theta_i^{-1}\)

Esta distribuição é imprópria porque a integral sobre o simplex \(\{\theta \in \mathbb{R}^k : \theta_i > 0, \sum_{i=1}^k \theta_i = 1\}\) não converge.

A função de verossimilhança da distribuição multinomial é:

\(p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}\)

Temos assim a posteriori:

\(p(\theta | n) \propto p(n | \theta) p(\theta) \propto \left( \prod_{i=1}^k \theta_i^{n_i} \right) \left( \prod_{i=1}^k \theta_i^{-1} \right) = \prod_{i=1}^k \theta_i^{n_i - 1}\)

Esta é a forma da distribuição Dirichlet com parâmetros \((n_1, n_2, \ldots, n_k)\):

\(p(\theta | n) \propto \prod_{i=1}^k \theta_i^{n_i - 1}\)

Para que a distribuição a posteriori seja própria, a integral da densidade da posteriori sobre o simplex deve ser finita. A distribuição Dirichlet \(\text{Dirichlet}(n_1, n_2, \ldots, n_k)\) é própria se e somente se todos os \(n_i > 0\). Se qualquer \(n_i = 0\), o termo \(\theta_i^{-1}\) aparece na densidade, e a integral não é finita.

Portanto, a distribuição a posteriori é própria se e somente se \(n_i > 0\) para todos os \(i = 1, 2, \ldots, k\).

Questão 3

Na véspera do primeiro turno para a eleição de governador do DF de 2010, a Datafolha divulgou uma pesquisa indicando que, de 891 eleitores entrevistados que já tinham decidido em quem votar, Agnelo Queiroz tinha a preferência de 467, Weslian Roriz a de 315 e outros candidatos a de 109 eleitores.

Formule um modelo para analisar esses dados. O interesse centra fundamentalmente em três perguntas:

  1. A eleição poderia ser definida no primeiro turno?

Usaremos um modelo baseado na distribuição Multinomial para representar os votos dos eleitores. A distribuição Multinomial é adequada aqui porque modela a probabilidade de um certo número de sucessos em várias categorias.

Seja \(\theta = (\theta_1, \theta_2, \theta_3)\) as proporções de votos para Agnelo, Weslian e outros candidatos, respectivamente, onde:

\(\theta_1 + \theta_2 + \theta_3 = 1\)

Os dados \(n = (n_1, n_2, n_3)\) seguem uma distribuição multinomial:

\((n_1, n_2, n_3) \sim \text{Multinomial}(n, \theta)\)

Para a análise bayesiana, assumimos uma distribuição a priori para \(\theta\). Podemos usar uma distribuição de Dirichlet com parâmetros \(\alpha = (\alpha_1, \alpha_2, \alpha_3)\):

\(\theta \sim \text{Dirichlet}(\alpha)\)

Dado que os dados seguem uma distribuição multinomial e a priori é uma Dirichlet, a distribuição a posteriori de \(\theta\) também é uma Dirichlet com parâmetros \(\alpha' = (\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)\):

\(\theta | n \sim \text{Dirichlet}(\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)\)

Para uma priori não informativa, podemos assumir \(\alpha_i = 1\) para todos os \(i\).

Para que a eleição seja definida no primeiro turno, um candidato deve obter mais de 50% dos votos válidos. Vamos calcular a probabilidade posteriori de \(\theta_1 > 0.5\):

\(P(\theta_1 > 0.5 | n)\)

  1. O candidato Agnelo poderia ser eleito no primeiro turno?

Semelhante à (a), precisamos calcular:

\(P(\theta_1 > 0.5 | n)\)

  1. Qual será a diferença na porcentagem dos votos válidos entre os dois primeiros colocados?

Calculamos a diferença na proporção dos votos válidos entre Agnelo e Weslian:

\(P(|\theta_1 - \theta_2| > d | n)\)

onde \(d\) é a diferença que desejamos calcular.

Podemos implementar isso usando simulações a partir da distribuição posteriori Dirichlet para estimar essas probabilidades.

Show/hide code
n_total = 891
n_agnelo = 467
n_weslian = 315
n_outros = 109

alpha = np.array([1, 1, 1])

alpha_post = alpha + np.array([n_agnelo, n_weslian, n_outros])

n_samples = 100000
samples = dirichlet(alpha_post).rvs(n_samples)

prob_agnelo_eleito = np.mean(samples[:, 0] > 0.5)
dif_agnelo_weslian = samples[:, 0] - samples[:, 1]
diff_percent = np.mean(dif_agnelo_weslian)

print(f"Probabilidade de Agnelo ser eleito no primeiro turno: {prob_agnelo_eleito:.2f}")
print(f"Diferença média na porcentagem entre Agnelo e Weslian: {diff_percent:.2f}")

plt.figure(figsize=(5,5))
x = np.linspace(0, 1, 1000)
plt.plot(x, beta.pdf(x, alpha_post[0], alpha_post[1] + alpha_post[2]), label="Agnelo")
plt.plot(x, beta.pdf(x, alpha_post[1], alpha_post[0] + alpha_post[2]), label="Weslian")
plt.axvline(0.5, color='red', linestyle='--', label="50% Threshold")
plt.xlabel("Proporção")
plt.ylabel("Densidade de Probabilidade")
plt.title("Distribuições Posteriori das Proporções de Votos")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper right')
plt.grid(True)
plt.tight_layout() 
plt.show()
Probabilidade de Agnelo ser eleito no primeiro turno: 0.92
Diferença média na porcentagem entre Agnelo e Weslian: 0.17

Com base nos dados da pesquisa e na análise bayesiana, é altamente provável que Agnelo Queiroz seja eleito no primeiro turno, com uma vantagem média de 17% sobre Weslian Roriz. A probabilidade de Agnelo obter mais de 50% dos votos válidos é muito alta (91.97%), o que sugere que a eleição pode ser decidida no primeiro turno.

Lista 5

Questão 2

Considere \(n = 12\) ensaios de Bernoulli com probabilidade de sucesso \(p\) e \(y = 9\) sucessos. Suponha que a distribuição a priori de \(p\) é especificada de forma que \(\eta = \log \left( \frac{p}{1 - p} \right)\) segue uma distribuição Normal com média \(\mu = 0\) e variância \(\sigma^2 = 1\).

De imediato temos:

\(\eta \sim N(0, \sigma^2)\)

O parâmetro \(p\) pode ser expresso em termos de \(\eta\) como:

\(p = \frac{e^\eta}{1 + e^\eta}\)

A verossimilhança para a distribuição Bernoulli é:

\(p(y | p) = \binom{n}{y} p^y (1 - p)^{n - y}\)

Para \(n = 12\) e \(y = 9\):

\(p(9 | p) = \binom{12}{9} p^9 (1 - p)^3\)

Para calcular \(\mathbb{E}(p | y = 9)\) e \(\text{DP}(p | y = 9)\), podemos usar a aproximação de Laplace, que envolve encontrar a moda da posteriori e aproximar a distribuição posteriori por uma Normal em torno da moda.

Podemos implementar isso em Python para calcular as aproximações solicitadas.

Show/hide code
n = 12
y = 9
mu = 0

variances = [1, 4, 9]

def posterior_approximation(y, n, mu, sigma2):

    # Log-verossimilhança em termos de eta
    def log_likelihood_eta(eta):
        p = expit(eta)  
        return y * np.log(p) + (n - y) * np.log(1 - p)
    
    def prior_eta(eta):
        return norm.logpdf(eta, loc=mu, scale=np.sqrt(sigma2))
    
    def log_posterior_eta(eta):
        return log_likelihood_eta(eta) + prior_eta(eta)
    
    def dlog_posterior_eta(eta):
        p = expit(eta)
        return y - n * p - eta / sigma2
    
    def d2log_posterior_eta(eta):
        p = expit(eta)
        return -n * p * (1 - p) - 1 / sigma2
    
    # Modo da posteriori (máximo a posteriori)
    from scipy.optimize import minimize
    result = minimize(lambda eta: -log_posterior_eta(eta), x0=0, method='Newton-CG', jac=lambda eta: -dlog_posterior_eta(eta), hess=lambda eta: -d2log_posterior_eta(eta))
    eta_mode = result.x[0]
    
    # Variância em torno da moda
    eta_var = -1 / d2log_posterior_eta(eta_mode)
    
    eta_samples = np.random.normal(loc=eta_mode, scale=np.sqrt(eta_var), size=10000)
    p_samples = expit(eta_samples)
    
    mean_p = np.mean(p_samples)
    std_p = np.std(p_samples)
    prob_p_greater_half = np.mean(p_samples > 0.5)
    
    return p_samples, mean_p, std_p, prob_p_greater_half

fig, axes = plt.subplots(len(variances), 1, figsize=(5,10))

for ax, sigma2 in zip(axes, variances):
    p_samples, mean_p, std_p, prob_p_greater_half = posterior_approximation(y, n, mu, sigma2)
    
    ax.hist(p_samples, bins=50, density=True, alpha=0.75)
    ax.axvline(mean_p, color='r', linestyle='--', linewidth=2, label=f'Média = {mean_p:.3f}')
    ax.axvline(0.5, color='g', linestyle='-.', linewidth=2, label='$p = 0.5$')
    ax.set_title(f'Posterior de $p$ para $\sigma^2 = {sigma2}$', fontsize=14)
    ax.set_xlabel('$p$', fontsize=12)
    ax.set_ylabel('Densidade', fontsize=12)
    ax.legend(loc='upper left', fontsize=12)
    ax.grid(True)
    
    
    print(f"Para σ² = {sigma2}:")
    print(f"  E(p | y = 9) = {mean_p:.4f}")
    print(f"  DP(p | y = 9) = {std_p:.4f}")
    print(f"  Pr(p > 0.5 | y = 9) = {prob_p_greater_half:.4f}")
    print()

plt.tight_layout()
plt.show()
Para σ² = 1:
  E(p | y = 9) = 0.6736
  DP(p | y = 9) = 0.1099
  Pr(p > 0.5 | y = 9) = 0.9275

Para σ² = 4:
  E(p | y = 9) = 0.7144
  DP(p | y = 9) = 0.1168
  Pr(p > 0.5 | y = 9) = 0.9502

Para σ² = 9:
  E(p | y = 9) = 0.7240
  DP(p | y = 9) = 0.1208
  Pr(p > 0.5 | y = 9) = 0.9469

Percebemos que a medida que a variância da priori aumenta, a incerteza na posteriori também aumenta, refletindo uma maior incerteza a priori sobre \(p\), sendo propagada para a distribuição a posteriori.